Usaremos una muestra aleatoria de tamaño 200 de la Encuesta nacional de Consumo Alimentario del 2014, ya trabajadas en clases anteriores (base_enca2014_200.xlsx)

Llamamos la librería readxl y cargamos nuestro excel

library(readxl)
datos <- read_excel("data/base_enca2014_200.xlsx")  #Cargamos la base

Agregamos las etiquetas a los códigos de las variables categóricas

datos$nse <- factor(datos$nse, labels = c("Alto", "Medio_alto", "Medio", "Medio_Bajo", "Bajo"))
datos$sex <- factor(datos$sex, labels = c("Hombre","Mujer"))
datos$macrozona <- factor(datos$macrozona, labels = c("Norte","Centro Norte","Centro Sur","Sur","RM"))
datos$area <- factor(datos$area, labels = c("Urbano","Rural"))

Seleccionamos la variable Peso en Kilogramos (wgt) y hacemos un summary() para revisar la variable

summary(datos$wgt)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  34.10   61.00   68.90   71.21   79.03  114.20 

Al no ver nada extraño, procedemos a verificar si esta variable si se ajusta a la normalidad

library(nortest)         #llamamos esta librería para usar el lillies.test()
lillie.test(datos$wgt)   #Aplicamos la prueba de normalidad de Kolmogorov-Smirnov (c/corrección de Lilliefors)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  datos$wgt
D = 0.076351, p-value = 0.006484

shapiro.test(datos$wgt)  #Aplicamos la prueba de normalidad de Shapiro-Wilk

    Shapiro-Wilk normality test

data:  datos$wgt
W = 0.97193, p-value = 0.0004919

Observamos que tanto el lillies.test como shapiro.test, rechazan la normalidad, por tanto procedemos a transformar la variable para evualuar nuevamente el ajuste.

#Aplicamos transformación de la raiz cuadrada
lillie.test(sqrt(datos$wgt))  #realizamos la transformación para garantizar el ajsute a la normal

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  sqrt(datos$wgt)
D = 0.057148, p-value = 0.1145

#shapiro.test(sqrt(datos$wgt)) #realizamos la transformación para garantizar el ajsute a la normal

Vemos que con la transformación de la raiz cuadrada, si se ajusta a la normal, por tanto, procedemos a realizar la transformación en la base para poder usarla posteriormente.

datos$wgt_sqrt <- sqrt(datos$wgt)

La nueva variable wgt_sqrt es la que vamos utilizar y procedemos a graficarlo

library(ggplot2)
ggplot(data=datos, aes(sample = wgt_sqrt)) + 
  stat_qq() + stat_qq_line(col="red")+
  xlab('Teoricas')+ylab('Observadas')+theme_light() + ggtitle("qq-plot de Peso transformado (wgt_sqrt)")

CONTRASTE DE HIPÓTESIS PARA LA MEDIA

Supongamos que conocemos que el promedio poblacional del peso adulto de los chilenos es 74.1 kgr producto del informe del Estado Nutricional del Ministerio de la Salud del 2017 y queremos saber si difiere del promedio muestral al 95% de confianza.

#Para aplicar el z-test es necesario instalar el paquete BSDA
#install.packages("BSDA")
#OJO: como la variable peso está transformada, el valor del parámetro poblacional 
#tambien debe estar transformado
library(BSDA)
z.test(x=datos$wgt_sqrt, mu=sqrt(74.1), sigma.x=sd(datos$wgt_sqrt), conf.level=0.95)

    One-sample z-Test

data:  datos$wgt_sqrt
z = -3.5601, p-value = 0.0003706
alternative hypothesis: true mean is not equal to 8.608136
95 percent confidence interval:
 8.280608 8.513191
sample estimates:
mean of x 
   8.3969 

En conclusión como el p-value es menor a 0.05, rechazamos la Ho de igualdad. Es decir, existen diferencias significativas entre el promedio muestral y el promedio poblacional.

Ahora aplicamos el t-test

t.test(x=datos$wgt_sqrt, mu=sqrt(74.1),conf.level = 0.95)

    One Sample t-test

data:  datos$wgt_sqrt
t = -3.5601, df = 199, p-value = 0.0004636
alternative hypothesis: true mean is not equal to 8.608136
95 percent confidence interval:
 8.279897 8.513903
sample estimates:
mean of x 
   8.3969 

Obtenemos la misma conclusión que con el z.test, es decir, rechazamos la Ho de igualdad, dado que el p-value es menor a 0.05, por tanto, existen diferencias significativas entre el promedio muestral y el promedio poblacional.

CONTRASTE PARA LA DIFERENCIA DE 2 MEDIAS MUESTRALES

Ahora queremos ver si hay diferencias en el peso por género. Hacemos un boxplot para visualizar.

ggplot(data=datos, aes(x=sex, y=wgt_sqrt, fill=sex)) +
  geom_boxplot(outlier.colour = "red" ) + 
  ggtitle('Grafico de Boxplot de Peso transformado por Genero') +
  stat_boxplot(geom = "errorbar", width = 0.25) + 
  stat_summary(fun=mean,geom="point",shape=18,size=3,color="black")+ theme_bw()

Para aplicar el z.test debemos segmentar el peso por el género

peso_mas <- datos[datos$sex=="Hombre" , "wgt_sqrt"]
peso_fem <- datos[datos$sex=="Mujer" , "wgt_sqrt"]

z.test(x=peso_mas$wgt_sqrt, sigma.x = sd(peso_mas$wgt_sqrt), 
       y=peso_fem$wgt_sqrt, sigma.y = sd(peso_fem$wgt_sqrt), conf.level=0.95)

    Two-sample z-Test

data:  peso_mas$wgt_sqrt and peso_fem$wgt_sqrt
z = 3.8338, p-value = 0.0001262
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2230608 0.6896987
sample estimates:
mean of x mean of y 
 8.704956  8.248577 

Conclusión: como el p-value es menor a 0.05, rechazamos la Ho de igualdad, es decir, el promedio del peso es diferente según el género. En otras palabras existe dimorfismo para el peso.

Para usar el t.test, debemos usar previamente la prueba de Levene para conocer si existe homocedasticidad entre los grupos. Para esto hay que instalar el paquete car.

library(car)
leveneTest(datos$wgt_sqrt, datos$sex) 

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   1  0.1844 0.6681
      198               

El test no permite rechazar la Ho, por tanto, decimos que existe homocedasticidad y especificamos var.equal=T para obtener un t de Student regular.

t.test(datos$wgt_sqrt~datos$sex,var.equal=TRUE ,conf.level=0.95)

    Two Sample t-test

data:  datos$wgt_sqrt by datos$sex
t = 3.7168, df = 198, p-value = 0.0002625
alternative hypothesis: true difference in means between group Hombre and group Mujer is not equal to 0
95 percent confidence interval:
 0.2142404 0.6985191
sample estimates:
mean in group Hombre  mean in group Mujer 
            8.704956             8.248577 

#t.test(peso_mas$wgt_sqrt, peso_fem$wgt_sqrt,var.equal=TRUE ,conf.level=0.95)

CONTRASTE DE HIPOTESIS PARA PROPORCIONES

Supongamos que conocemos el porcentaje poblacional para el consumo Energía proveniente de carbohidratos en 51%

Para obtener el porcentaje de la muestra usaremos la variable ChoEnerPor. Como cada valor es un porcentaje, usaremos el promedio geométrico, ya previamente explicado, para obtener un único porcentaje.

library(psych)  # para usar el geometric.mean(), hay que usar este paquete
geometric.mean(datos$ChoEnerPor)

[1] 58.79112

El porcentaje de 58.79112, tenemos que convertirlo a proporción para poder convertirlo en un frecuencia absoluta.

#La frecuencia la guardamos en un objeto llamado cho
cho <- round((geometric.mean(datos$ChoEnerPor) / 100)*nrow(datos))
cho

[1] 118

La frecuencia absoluta es 118, es decir hay 118 personas que consumen carbohidratos.

CONTRASTE DE HIPOTESIS PARA PROPORCIONES

Para realizar el contraste de proporciones usaremos la función prop.test(), previamente ya descrita en clases anteriores.

#OJO: El parámetro porblacional se debe transformar a proporción.
prop.test(x=cho, n=nrow(datos), p=0.51, conf.level=0.95)  

    1-sample proportions test with continuity correction

data:  cho out of nrow(datos), null probability 0.51
X-squared = 4.8069, df = 1, p-value = 0.02835
alternative hypothesis: true p is not equal to 0.51
95 percent confidence interval:
 0.5182510 0.6582337
sample estimates:
   p 
0.59 

Conclusión: como el p-value es menor a 0.05, rechazamos la Ho de igualdad, esto quiere decir que existe diferencia significativa de la proporción entre la muestra y la población

CONTRASTE DE HIPOTESIS PARA DIFERENCIA DE PROPORCIONES

Ahora queremos saber si existe diferencia por género entre la porporción de los consumidores de carbohidratos.

Para esto, hacemos el mismo procedimiento, pero segmentandolo por género. Primero obtenemos los “n” de cada género

n_mas <- nrow(datos[datos$sex=="Hombre" , "ChoEnerPor"])  #65
n_fem <- nrow(datos[datos$sex=="Mujer" , "ChoEnerPor"])  #135
n_mas

[1] 65

n_fem

[1] 135

Luego obtenemos las frecuencias absolutas del consumo de carbihidratos para cada género.

cho_mas <- round((geometric.mean(datos[datos$sex=="Hombre" , "ChoEnerPor"]) / 100)*n_mas) #hay 39
cho_fem <- round((geometric.mean(datos[datos$sex=="Mujer" , "ChoEnerPor"]) / 100)*n_fem) #hay 79
cho_mas

ChoEnerPor 
        39 

cho_fem

ChoEnerPor 
        79 

CONTRASTE DE HIPOTESIS PARA DIFERENCIA DE PROPORCIONES

Luego hacemos el contraste con los valores previamente calculados

prop.test(x=c(cho_mas, cho_fem),  n=c(n_mas,n_fem), conf.level=0.95)

    2-sample test for equality of proportions with continuity correction

data:  c(cho_mas, cho_fem) out of c(n_mas, n_fem)
X-squared = 0.00212, df = 1, p-value = 0.9633
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.1418092  0.1714388
sample estimates:
   prop 1    prop 2 
0.6000000 0.5851852 

Conclusión: como el p-value es mayor a 0.05, no podemos rechazar la Ho, por tanto decimos que no hay diferencias entre la porporción de consumidores de carbohidratos por género.

CORRELACION

Para aplicar la correlación de Pearson, debemos verificar la normalidad de las dos variables continuas a estudiar. Por ejemplo el Peso (wgt) y la talla (wgt)

lillie.test(sqrt(datos$wgt))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  sqrt(datos$wgt)
D = 0.057148, p-value = 0.1145

lillie.test(sqrt(datos$ht))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  sqrt(datos$ht)
D = 0.092522, p-value = 0.0002663

Como ninguna se ajusta a la normal, intentamos hacer transformaciones para mejorarlas.

shapiro.test(sqrt(datos$wgt))

    Shapiro-Wilk normality test

data:  sqrt(datos$wgt)
W = 0.98691, p-value = 0.06177

shapiro.test(sqrt(datos$ht))

    Shapiro-Wilk normality test

data:  sqrt(datos$ht)
W = 0.97031, p-value = 0.0003085

Como solo el peso se ajusta la normal, no podemos usar Pearson. Por tanto, podemos usar la correlación no-paramétrica: Spearman o kendall

CORRELACION

Comenzamos usando el coeficiente de correlación de Spearman

cor.test(datos$wgt, datos$ht, method = "spearman",exact=FALSE)

    Spearman's rank correlation rho

data:  datos$wgt and datos$ht
S = 739489, p-value = 3.893e-11
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.4453693 

Probamos tambien el coeficiente de correlación de kendall

cor.test(datos$wgt, datos$ht, method = "kendall",exact=FALSE)

    Kendall's rank correlation tau

data:  datos$wgt and datos$ht
z = 6.353, p-value = 2.111e-10
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.3035102 

CORRELACION

Graficamos la correlación, mediante el Grafico de Dispersión

ggplot(datos, aes(x = ht, y = wgt)) + 
  geom_point() + ggtitle('Grafico de Dispersión') + theme(legend.position = "none") +
  theme_light() + geom_smooth(method = "lm", se = FALSE)

MATRIZ DE CORRELACIONES

Buscar la correlación entre 3 variables wgt, energia_kc, cho_g_dia, Pero primero evaluamos si se ajustan a la normal.

lillie.test(log(datos$wgt))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  log(datos$wgt)
D = 0.037743, p-value = 0.6955

lillie.test(log(datos$energia_kc))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  log(datos$energia_kc)
D = 0.061803, p-value = 0.0607

lillie.test(log(datos$cho_g_dia))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  log(datos$cho_g_dia)
D = 0.048139, p-value = 0.3109

Hacemos una nueva base solo con las variables seleccionas

datos2 <- datos[, c("wgt", "energia_kc", "cho_g_dia")]
cor(log(datos2), method = "pearson")  #la función cor(), nos permite una matriz de correlación

                    wgt  energia_kc   cho_g_dia
wgt         1.000000000 0.001472497 -0.05981235
energia_kc  0.001472497 1.000000000  0.89986169
cho_g_dia  -0.059812347 0.899861691  1.00000000

Podemos hacer una representación gráfica de la correlación con el paquete corrplot

#install.packages("corrplot")    #Lo instalamos previamente
library(corrplot)

MATRIZ DE CORRELACIONES

Entonces con la función corrplot() podemos representar gráficamente una matriz de correlaciones

corrplot(cor(log(datos2), method = "pearson"), method ="number", 
         tl.col = "black", tl.srt = 20,
         col=colorRampPalette(c("blue","lightblue","red"))(100))

CORRELACIONES PARCIALES

Para realizar las correlaciones parciales, vamos a plantearnos la correlación entre el Peso (wgt) y la Talla (ht), controlando la Edad (ageyrs)

Sabemos que la talla no se ajusta a normal. Aplicaremos la correlación de Spearman

#Volvemos a realizar la correlación de spearman sin controlar el efecto de la edad
cor.test(datos$wgt, datos$ht, method = "spearman",exact=FALSE)  

    Spearman's rank correlation rho

data:  datos$wgt and datos$ht
S = 739489, p-value = 3.893e-11
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.4453693 

#install.packages("ppcor")    #Usamos paquete ppcor para realizar la correlación parcial
library(ppcor)
pcor.test(x=datos$ht, y=datos$wgt, z=datos$ageyrs, method = "spearman")  #Aqui controlamos por la edad

   estimate      p.value statistic   n gp   Method
1 0.4573582 1.113879e-11  7.218548 200  1 spearman

Observamos que la correlación aumentó un poco (0.457) y es significativa (p-value < 0.05)

Por tanto, la correlación aumenta si la edad permanece constante.